
# Import Libraries

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns
import plotly
import plotly.offline as pyoff
import plotly.graph_objs as go
import plotly.express as px
import chart_studio
import chart_studio.plotly as py
import calmap
import datetime
import tensorflow as tf
import os
import random
import re
import plotly.offline as pyoff
import plotly.graph_objs as go

from datetime import date
from plotly.subplots import make_subplots
from itertools import cycle, product
from statsmodels.tsa.seasonal import STL
from scipy.stats import boxcox
from pmdarima.arima import ADFTest, KPSSTest, auto_arima
from pmdarima.utils import diff_inv
from statsmodels.tsa.stattools import adfuller
from sklearn.model_selection import TimeSeriesSplit
from tensorflow.keras.layers import LSTM, Dense, BatchNormalization
from tensorflow.keras import Sequential
from tensorflow.keras.backend import clear_session
from tensorflow.keras.callbacks import EarlyStopping
from tensorflow.keras.preprocessing.sequence import TimeseriesGenerator
from tensorflow.keras.initializers import *
from tensorflow.keras import optimizers
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from sklearn.linear_model import LinearRegression
from scipy.special import boxcox1p, inv_boxcox1p
import matplotlib.patches as mpatches
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from sklearn.model_selection import GridSearchCV
from joblib import delayed
from warnings import catch_warnings
from warnings import filterwarnings
from statsmodels.tsa.forecasting.stl import STLForecast
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.preprocessing import StandardScaler
from tensorflow.keras.optimizers import Adam
from IPython.core.display import HTML, display


# Versões dos pacotes usados neste jupyter notebook
%reload_ext watermark
%watermark -a "Herikc Brecher" --iversions

Author: Herikc Brecher

plotly      : 4.14.3
re          : 2.2.1
tensorflow  : 2.5.0
seaborn     : 0.11.1
pandas      : 1.2.2
chart_studio: 1.1.0
datetime    : 4.3
numpy       : 1.19.5
calmap      : 0.0.9
matplotlib  : 3.4.2


# Variaveis globais
SEED = 84796315
PREVISOES = 5
EPOCAS = 30
BATCH_SIZE = 1000
EXECUTAR_GRID_SEARCH = False


# Configurando seeds
os.environ['PYTHONHASHSEED'] = str(SEED)
tf.random.set_seed(SEED)
np.random.seed(SEED)
random.seed(SEED)


# Exibindo toda tela
display(HTML("<style>.container { width:100% !important; }</style>"))


# Import dataset
dtOrders = pd.read_csv('../data/olist_orders_dataset.csv', encoding = 'utf8')


# Colunas do tipo data
dateColumns = ['order_purchase_timestamp', 'order_approved_at', 'order_delivered_carrier_date',\
               'order_delivered_customer_date', 'order_estimated_delivery_date']

# Dataset de analise temporal
dtOrdersAjustado = dtOrders.copy()


# Convertendo columas de data para date
for col in dateColumns:
    dtOrdersAjustado[col] = pd.to_datetime(dtOrdersAjustado[col], format = '%Y-%m-%d %H:%M:%S')


# Dropando valores NA
dtOrdersAjustado = dtOrdersAjustado.dropna()


dtOrdersAjustado.dtypes

order_id                                 object
customer_id                              object
order_status                             object
order_purchase_timestamp         datetime64[ns]
order_approved_at                datetime64[ns]
order_delivered_carrier_date     datetime64[ns]
order_delivered_customer_date    datetime64[ns]
order_estimated_delivery_date    datetime64[ns]
dtype: object


dtHistory = pd.to_datetime(dtOrdersAjustado['order_purchase_timestamp']).dt.date

start = dtHistory.min()
end = dtHistory.max()

idx = pd.date_range(start, end, normalize=True)

seriesHistory = dtHistory.value_counts(sort=False).sort_index().reindex(idx, fill_value=0)

dtHistory = pd.DataFrame(seriesHistory).reset_index()


dtHistory


# Plot

# Definição dos dados no plot (Iniciando em Fevereiro de 2017 para não destorcer os dados)
plot_data = [go.Scatter(x = dtHistory['index'],
                        y = dtHistory['order_purchase_timestamp'])]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Periodo'},
                        yaxis = {"title": 'Vendas'},
                        title = 'Vendas por dia')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

pyoff.iplot(fig)


# Remove outliers
seriesHistory = seriesHistory[datetime.date(2017, 1, 1): datetime.date(2018, 8, 17)]
pred_range = pd.date_range(datetime.date(2018, 8, 17), datetime.date(2018, 10, 17))
dtHistory = pd.DataFrame(seriesHistory).reset_index()


# Plot

# Definição dos dados no plot (Iniciando em Fevereiro de 2017 para não destorcer os dados)
plot_data = [go.Scatter(x = dtHistory['index'],
                        y = dtHistory['order_purchase_timestamp'])]

# Layout
plot_layout = go.Layout(xaxis = {'title': 'Periodo'},
                        yaxis = {"title": 'Vendas'},
                        title = 'Vendas por dia')

# Plot da figura
fig = go.Figure(data = plot_data, layout = plot_layout)

pyoff.iplot(fig)


#Plot histórico de vendas por dia, mês e ano
fig, caxs = calmap.calendarplot(seriesHistory, daylabels='MTWTFSS', fillcolor='grey',cmap='YlGn', fig_kws=dict(figsize=(18, 9)))
fig.suptitle('Histórico de Vendas', fontsize=22)

fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.03, 0.67])
fig.colorbar(caxs[0].get_children()[1], cax=cbar_ax)

plt.show()


# Criar grafico na estrutura STL 4 layers
def add_stl_plot(fig, res, legend):
    axs = fig.get_axes()
    
    # Nome de cada um dos subplots
    comps = ["trend", "seasonal", "resid"]
    for ax, comp in zip(axs[1:], comps):
        series = getattr(res, comp)
        if comp == "resid":
            ax.plot(series, marker="o", linestyle="none")
        else:
            ax.plot(series)
            ax.legend(legend, frameon=False)


# Gerar STL
stl = STL(seriesHistory)
stl_res = stl.fit()
fig = stl_res.plot()
fig.set_size_inches((20, 12))
plt.show()


# Gerar STL não robusto e concatenar ao robusto
stl = STL(seriesHistory, robust=True)
res_robust = stl.fit()
fig = res_robust.plot()
fig.set_size_inches((20, 12))
res_non_robust = STL(seriesHistory, robust=False).fit()
add_stl_plot(fig, res_non_robust, ["Robusto", "Não Robusto"])


# Gerando STL para separar cada um dos componentes
stl = STL(seriesHistory)
res = stl.fit()

# Separando deseasonal
deseasonal = res.observed - res.seasonal

# Separando boxcox
bc_history, lmbda = boxcox(seriesHistory+1)
bc_history = pd.Series(bc_history, index=seriesHistory.index)

# Separando stationary
diff_history = seriesHistory.diff(7).dropna()
xi = seriesHistory.iloc[:7]

# Plotando todos elementos de forma separada
fig, axs = plt.subplots(nrows=4, sharex=True, figsize=(14, 8))
seriesHistory.plot(ax=axs[0])
axs[0].set_ylabel('Original')
deseasonal.plot(ax=axs[1])
axs[1].set_ylabel('Deseasonal')
bc_history.plot(ax=axs[2])
axs[2].set_ylabel('BoxCox')
diff_history.plot(ax=axs[3])
axs[3].set_ylabel('Stationary')
fig.align_ylabels()
fig.suptitle('Transformações')
plt.tight_layout()
plt.show()


# Teste ADF automatizado
def adf_test(series):
    result = adfuller(series, autolag='AIC')

    print(f'ADF: {result[0]}')
    print(f'steps: {result[1]}')
    print(f'p-value: {result[1]}')

    for key, value in result[4].items():
        print('\nValor Critico:')
        print(f'   {key}: {value}')


adf = ADFTest(alpha = 0.05)
adf.should_diff(seriesHistory)

(0.01, False)


# Testando estacionaridade na série original
adf_test(seriesHistory)

ADF: -2.6163082378564684
steps: 0.08967331201291256
p-value: 0.08967331201291256

Valor Critico:
   1%: -3.441694608475642

Valor Critico:
   5%: -2.866544718556839

Valor Critico:
   10%: -2.5694353738653684


adf.should_diff(diff_history)

(0.01, False)


adf_test(diff_history)

ADF: -6.801606152090551
steps: 2.229745766641977e-09
p-value: 2.229745766641977e-09

Valor Critico:
   1%: -3.441834071558759

Valor Critico:
   5%: -2.8666061267054626

Valor Critico:
   10%: -2.569468095872659


kpss = KPSSTest(alpha = 0.05)

# Testando estacionaridade na série original com KPSS
kpss.should_diff(seriesHistory)

(0.01, True)


kpss.should_diff(diff_history)

(0.1, False)


# Controle de resultados de toda fase de modelagem
result = pd.DataFrame(columns=['Algorithm', 'MSE', 'RMSE', 'MAE', 'Mean_Real_Value', 'Mean_Predict_Value'])


split_range = TimeSeriesSplit(n_splits = 8, max_train_size = pred_range.shape[0], test_size = pred_range.shape[0])


# Adiciona o registro ao dataset
def record(result, algorithm, mse = -1, rmse = -1, mae = -1, mrv = -1, mpv = -1, show = True):
    new = pd.DataFrame(dict(Algorithm = algorithm, MSE = mse, RMSE = rmse, MAE = mae, Mean_Real_Value = mrv,\
                            Mean_Predict_Value = mpv), index=[0])
    result = pd.concat([result, new], ignore_index=True)
    
    if show:
        display(result)
    
    return result


# Plot no formato de 4 layers, seguindo o STL para cada um dos modelos
def plot(index, pred, mse, title, fig = None, ax = None, ylabel = ''):
    global seriesHistory
    
    empty_fig = fig is None
    
    if empty_fig:
        fig, ax = plt.subplots(figsize=(13, 6))
    else: 
        ax.set_ylabel(ylabel)
                
    ax.set_title(title)    
    patch_ = mpatches.Patch(color = 'white', label = f'MSE: {np.mean(mse):.1e}')
    L1 = ax.legend(handles = [patch_], loc = 'upper left', fancybox = True, framealpha = 0.7,  handlelength = 0)
    ax.add_artist(L1)
    
    sns.lineplot(x = seriesHistory.index, y = seriesHistory, label = 'Real', ax = ax)
    sns.lineplot(x = index, y = pred, label = 'Previsto', ax = ax)
    ax.axvline(x = index[0], color = 'red')
    ax.legend(loc = 'upper right')
    
    if empty_fig:
        plt.show()
    else:
        return fig


'''
Correção da função diff_inv original que apresenta problemas ainda não corrigidos:
https://github.com/alkaline-ml/pmdarima/issues/410
'''

def diff_inv_fix(x_diff, xi, column, steps = 1):
    total_len = len(x_diff) + len(xi)
    
    ix = pd.date_range(xi.index[0], periods = total_len)  
    inv = diff_inv(x_diff, steps, xi = xi) + np.fromiter(cycle(xi), count = total_len, dtype = float)  
    inv = pd.Series(inv, index = ix, name = column)
    
    return inv


# Calculo para previsão e teste quando utilizado a série Original
def calcPredTestOriginal(train, pred, test):
    return pred, test, 0


# Calculo para previsão e teste quando utilizado a série Deseasonal
def calcPredTestDeseasonal(train, pred, test):
    # Removendo a sazonalidade da série e convertendo para o shape correto
    last_seasonal = res.seasonal.reindex_like(train).tail(stl.period)
    pred = pred + np.fromiter(cycle(last_seasonal), count = pred.shape[0], dtype = float)
    test = test + res.seasonal.reindex_like(test)
    
    return pred, test, 1


# Calculo para previsão e teste quando utilizado a série BoxCox
def calcPredTestBoxCox(train, pred, test):
    # Reverdendo a normalização do boxcox
    pred = inv_boxcox1p(pred, lmbda)
    test = inv_boxcox1p(test, lmbda)
    
    return pred, test, 2


# Calculo para previsão e teste quando utilizado a série Stationary
def calcPredTestStationary(train, pred, test):
    # Calculando a diferença da sazonalidade
    xi = seriesHistory.reindex_like(train).tail(PREVISOES)
    pred = diff_inv_fix(pred, xi, 'order_purchase_timestamp', PREVISOES).iloc[PREVISOES:]
    test = diff_inv_fix(test, xi, 'order_purchase_timestamp', PREVISOES).iloc[PREVISOES:]
    
    return pred, test, 3


# Report para Time Series Regressor, realiza o treino do modelo, adiciona aos resultados e faz o plot de acompanhamento
def reportTSR(data, modelName, calcFunction):
    global result
    global figs
    
    mse = []
    rmse = []
    mae = []
    mrv = []
    mpv = []
    
    title = modelName + ' - Time Series Regression'
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
    
        gen = TimeseriesGenerator(train, train, PREVISOES, batch_size = BATCH_SIZE)

        X_train = gen[0][0]
        y_train = gen[0][1]

        lr = LinearRegression()
        lr.fit(X_train, y_train)
        X_pred = y_train[-PREVISOES:].reshape(1,-1)
        pred = np.empty(test.shape[0])

        for i in range(len(pred)):
            forecast = lr.predict(X_pred)
            X_pred = np.delete(X_pred, 0, 1)
            X_pred = np.concatenate((X_pred, forecast.reshape(-1, 1)), 1)    
            pred[i] = forecast
        
        pred, test, indexPlot = calcFunction(train, pred, test)

        mse.append(mean_squared_error(pred, test, squared = True))
        rmse.append(mean_squared_error(pred, test, squared = False))
        mae.append(mean_absolute_error(pred, test))
        mrv.append(np.mean(test))
        mpv.append(np.mean(pred))
    
    result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mrv), np.mean(mpv), False)
    return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)


# Reset da figura
figs, axs = plt.subplots(nrows = 4, sharex = True, figsize = (13,6))
figs.tight_layout()
plt.close()


reportTSR(seriesHistory.copy(), 'Original', calcPredTestOriginal)


reportTSR(deseasonal.copy(), 'Deseasonal', calcPredTestDeseasonal)


reportTSR(bc_history.copy(), 'BoxCox', calcPredTestBoxCox)


reportTSR(diff_history.copy(), 'Stationary', calcPredTestStationary)


result


# Função utilizada para o hypertuning de alpha, beta e gamma do Exponential Smoothing
def GSES(data, modelName, alpha, beta, gamma, calcFunction):    
    mse = []
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
        
        try:
            with catch_warnings():
                filterwarnings("ignore")
                ES = (
                    ExponentialSmoothing(train, trend='add', seasonal='add', seasonal_periods=PREVISOES)
                    .fit(smoothing_level = alpha, smoothing_trend = beta, smoothing_seasonal = gamma, method='ls')
                )

                pred = ES.forecast(test.shape[0])

                pred, test, _ = calcFunction(train, pred, test)

                mse.append(mean_squared_error(pred, test, squared = True))
        
        except:
            mse.append(-1)
    
    return np.mean(mse)


# Função utilizada para o hypertuning de demais parâmetros do Exponential Smoothing
def GSESOPT(data, modelName, trend, season, periods, bias, method, calcFunction):
    mse = []
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
        
        try:
            with catch_warnings():
                filterwarnings("ignore")
                ES = (
                    ExponentialSmoothing(train, trend = trend, seasonal = season, seasonal_periods = periods)
                    .fit(remove_bias = bias, method = method, optimized = True)
                )

                pred = ES.forecast(test.shape[0])

                pred, test, _ = calcFunction(train, pred, test)

                mse.append(mean_squared_error(pred, test, squared = True))       
        except:
            mse.append(-1)
    
    return np.mean(mse)


# Report para Exponential Smoothing, realiza o treino do modelo, adiciona aos resultados e faz o plot de acompanhamento
def reportES(data, modelName, model_kwargs, fit_kwargs, calcFunction):
    global result
    global figs
    
    mse = []
    rmse = []
    mae = []
    mrv = []
    mpv = []
    
    title = modelName + ' - Exponential Smoothing'
    indexPlot = 0
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
        
        ES = ExponentialSmoothing       
        stlf = STLForecast(train, ES, model_kwargs = model_kwargs).fit(fit_kwargs = fit_kwargs)
        pred = stlf.forecast(test.shape[0])
             
        pred, test, indexPlot = calcFunction(train, pred, test)

        mse.append(mean_squared_error(pred, test, squared = True))
        rmse.append(mean_squared_error(pred, test, squared = False))
        mae.append(mean_absolute_error(pred, test))
        mrv.append(np.mean(test))
        mpv.append(np.mean(pred))
    
    result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mrv), np.mean(mpv), False)
    return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)


# Função para gerar tabela de hypertuning ampla
def exp_smoothing_configs(seasonal=[None]):
    models = list()
    # Lista de argumentos
    t_params = ['add', 'mul', None]
    s_params = ['add', 'mul', None]
    p_params = seasonal
    r_params = [True, False]
    method_params = ['L-BFGS-B' , 'TNC', 'SLSQP', 'Powell', 'trust-constr', 'bh', 'ls']
    
    # Gerando lista de argumentos
    for t in t_params:
        for s in s_params:
            for p in p_params:
                for r in r_params:
                    for m in method_params:
                        cfg = [t, s, p, r, m]
                        models.append(cfg)
    return models


# Gerando tabela de hypertunning
alphas = betas = gammas = np.arange(1, step=0.1)
tuning = pd.DataFrame(product(alphas, betas, gammas), columns=['alpha', 'beta', 'gamma'])
tuning.head()


%%time

# Treinamento do modelo 
tuning['mse'] = tuning.apply(lambda x: GSES(seriesHistory.copy(), 'Original',\
                                            x.alpha, x.beta, x.gamma, calcPredTestOriginal), axis = 1)

Wall time: 55.5 s


# Verificando o menor mse
tuning.query('mse == mse.min() and mse != -1')


# Criando lista de argumentos ampla
params_ = exp_smoothing_configs([PREVISOES])


tuning2 = pd.DataFrame(params_, columns=['trend', 'season', 'periods', 'bias', 'method'])


len(tuning2)

126


tuning2.head()


%%time

# Se True irá treinar com a nova lista mais ampla (pode demorar)
if EXECUTAR_GRID_SEARCH:
    tuning2['mse'] = tuning2.apply(lambda x: GSESOPT(seriesHistory.copy(), 'Original',\
                                                     x.trend, x.season, x.periods, x.bias, x.method, calcPredTestOriginal),\
                                   axis = 1)

Wall time: 0 ns


if EXECUTAR_GRID_SEARCH:
    display(tuning2.query('mse == mse.min() and mse != -1'))


# Reset da figura
figs, axs = plt.subplots(nrows=4, sharex=True, figsize=(13,6))
figs.align_ylabels()
figs.tight_layout()
plt.close()


# Dicionar de configurações de acordo com o hypertuning
model_kwargs = dict(trend = None, seasonal = None, seasonal_periods = PREVISOES)
fit_kwargs = dict(remove_bias = False, smoothing_level = 0.1, smoothing_trend = 0.1, smoothing_seasonal = 0.3,\
                  method = 'Powell')


reportES(seriesHistory.copy(), 'Original', model_kwargs, fit_kwargs, calcPredTestOriginal)


alphas = betas = gammas = np.arange(1, step=0.1)
tuning = pd.DataFrame(product(alphas, betas, gammas), columns=['alpha', 'beta', 'gamma'])
tuning.head()


%%time
tuning['mse'] = tuning.apply(lambda x: GSES(deseasonal.copy(), 'Deseasonal',\
                                            x.alpha, x.beta, x.gamma, calcPredTestDeseasonal), axis=1)

Wall time: 58.4 s


tuning.query('mse == mse.min() and mse != -1')


params_ = exp_smoothing_configs([PREVISOES])


tuning2 = pd.DataFrame(params_, columns=['trend', 'season', 'periods', 'bias', 'method'])


%%time
if EXECUTAR_GRID_SEARCH:
    tuning2['mse'] = tuning2.apply(lambda x: GSESOPT(deseasonal.copy(), 'Deseasonal',\
                                                     x.trend, x.season, x.periods, x.bias, x.method, calcPredTestDeseasonal),\
                                   axis = 1)

Wall time: 0 ns


if EXECUTAR_GRID_SEARCH:
    display(tuning2.query('mse == mse.min() and mse != -1'))


model_kwargs = dict(trend = None, seasonal = 'add', seasonal_periods = PREVISOES)
fit_kwargs = dict(remove_bias = True, smoothing_level = 0.1, smoothing_trend = 0.1, smoothing_seasonal = 0.9, method = 'Powell')


reportES(deseasonal.copy(), 'Deseasonal', model_kwargs, fit_kwargs, calcPredTestDeseasonal)


alphas = betas = gammas = np.arange(1, step=0.1)
tuning = pd.DataFrame(product(alphas, betas, gammas), columns=['alpha', 'beta', 'gamma'])
tuning.head()


%%time
tuning['mse'] = tuning.apply(lambda x: GSES(bc_history.copy(), 'BoxCox',\
                                            x.alpha, x.beta, x.gamma, calcPredTestBoxCox), axis=1)

Wall time: 56.2 s


tuning.query('mse == mse.min() and mse != -1')


params_ = exp_smoothing_configs([PREVISOES])


tuning2 = pd.DataFrame(params_, columns=['trend', 'season', 'periods', 'bias', 'method'])


%%time
if EXECUTAR_GRID_SEARCH:
    tuning2['mse'] = tuning2.apply(lambda x: GSESOPT(bc_history.copy(), 'BoxCox',\
                                                  x.trend, x.season, x.periods, x.bias, x.method, calcPredTestBoxCox),\
                                   axis = 1)

Wall time: 0 ns


if EXECUTAR_GRID_SEARCH:
    display(tuning2.query('mse == mse.min() and mse != -1'))


model_kwargs = dict(trend='add', seasonal='add', seasonal_periods=PREVISOES)
fit_kwargs = dict(remove_bias = True, smoothing_level = 0.1, smoothing_trend = 0.1, smoothing_seasonal = 0.4,\
                  method = 'Powell')


reportES(bc_history.copy(), 'BoxCox', model_kwargs, fit_kwargs, calcPredTestBoxCox)


alphas = betas = gammas = np.arange(1, step=0.1)
tuning = pd.DataFrame(product(alphas, betas, gammas), columns=['alpha', 'beta', 'gamma'])
tuning.head()


%%time
tuning['mse'] = tuning.apply(lambda x: GSES(diff_history.copy(), 'Stationary',\
                                            x.alpha, x.beta, x.gamma, calcPredTestStationary), axis=1)

Wall time: 1min 2s


tuning.query('mse == mse.min() and mse != -1')


params_ = exp_smoothing_configs([PREVISOES])


tuning2 = pd.DataFrame(params_, columns=['trend', 'season', 'periods', 'bias', 'method'])


%%time
if EXECUTAR_GRID_SEARCH:
    tuning2['mse'] = tuning2.apply(lambda x: GSESOPT(diff_history.copy(), 'Stationary',\
                                                  x.trend, x.season, x.periods, x.bias, x.method, calcPredTestStationary),\
                                   axis = 1)

Wall time: 0 ns


if EXECUTAR_GRID_SEARCH:
    display(tuning2.query('mse == mse.min() and mse != -1'))


model_kwargs = dict(trend='add', seasonal='add', seasonal_periods=PREVISOES)
fit_kwargs = dict(remove_bias = False, smoothing_level = 0, smoothing_trend = 0.3, smoothing_seasonal = 0.0,\
                  method = 'ls')


reportES(diff_history.copy(), 'Stationary', model_kwargs, fit_kwargs, calcPredTestStationary)


result


# Report do algoritmo arima, também é adicionado a base de resultados e realizado o plot de acompanhamento
def reportArima(arimaModel, modelName, calcFunction):
    global result
    global figs
    
    mse = []
    rmse = []
    mae = []
    mrv = []
    mpv = []
    
    title = modelName + ' - '  + arimaModel.__str__().strip()
    indexPlot = 0
    
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
        arimaModel.fit(train)
        pred = arimaModel.predict(test.shape[0])
             
        pred, test, indexPlot = calcFunction(train, pred, test)

        mse.append(mean_squared_error(pred, test, squared = True))
        rmse.append(mean_squared_error(pred, test, squared = False))
        mae.append(mean_absolute_error(pred, test))
        mrv.append(np.mean(test))
        mpv.append(np.mean(pred))
    
    result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mrv), np.mean(mpv), False)
    return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)


# Reset da figura
figs, axs = plt.subplots(nrows=4, sharex=True, figsize=(13,6))
figs.align_ylabels()
figs.tight_layout()
plt.close()


# Correlação entre os periodos com ARIMA

lags = 90

with catch_warnings():
    filterwarnings("ignore")
    fig, ax = plt.subplots(2, figsize=(12, 6), sharex=True)
    plot_acf(seriesHistory.diff().dropna(), ax = ax[0], lags = lags, missing = 'drop')
    plot_pacf(seriesHistory.diff().dropna(), ax = ax[1], lags = lags)
    plt.show()


%%time

# Utilizando o auto arima para descobrir os argumentos ideias baseados no conjunto de dado informado
data = seriesHistory.copy()
arimaModel = auto_arima(seriesHistory.copy(), m = PREVISOES)
arimaModel

Wall time: 17.7 s

ARIMA(order=(2, 1, 1), scoring_args={}, seasonal_order=(1, 0, 1, 5),
      suppress_warnings=True, with_intercept=False)


reportArima(arimaModel, 'Original', calcPredTestOriginal)


%%time
data = deseasonal.copy()
arimaModel = auto_arima(data, seasonal = False)
arimaModel

Wall time: 2.66 s

ARIMA(order=(2, 1, 1), scoring_args={}, suppress_warnings=True)


reportArima(arimaModel, 'Deseasonal', calcPredTestDeseasonal)


%%time
data = bc_history.copy()
arimaModel = auto_arima(data, m = PREVISOES)
arimaModel

Wall time: 36.9 s

ARIMA(order=(4, 1, 2), scoring_args={}, seasonal_order=(0, 0, 0, 5),
      suppress_warnings=True)


reportArima(arimaModel, 'BoxCox', calcPredTestBoxCox)


%%time
data = diff_history.copy()
arimaModel = auto_arima(data, seasonal = False)
arimaModel

Wall time: 12.3 s

ARIMA(order=(3, 0, 3), scoring_args={}, suppress_warnings=True)


reportArima(arimaModel, 'Stationary', calcPredTestStationary)


result


# Redefinindo variaveis globais para o treino utilizando LSTM

BATCH_SIZE = 30
PREVISOES = 5


# Hypertuning do algoritmo de LSTM
def GSLSTM(data, learning_rate, calcFunction):
    mse = []
    
    # Crossvalidation para cada parte do conjunto
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]

        try:
            with catch_warnings():
                filterwarnings("ignore")
                
                # Normalização e reshape do conjunto de treino
                ss = StandardScaler()
                ss.fit(train.values.reshape(-1, 1))
                train_input = ss.transform(train.values.reshape(-1, 1))
                
                # Gerando conjunto de treino com TimeseriesGenerator baseado no conjunto atual e o batch informado
                test_input = train_input[-(PREVISOES + 1):]
                test_gen = TimeseriesGenerator(test_input, test_input, length = PREVISOES, batch_size = BATCH_SIZE)
                train_gen = TimeseriesGenerator(train_input, train_input, length = PREVISOES, batch_size = BATCH_SIZE)
                
                # Reset da sessão
                clear_session()
                
                # Construindo o modelo de LSTM com GlorotUniform pois inicializa de forma normalizada
                initializer = GlorotUniform(seed = SEED)
                model = Sequential()
                
                # 1 camada de LSTM com 64 entradas, 2 camadas densas e uma de normalização intermediando as camadas densas
                model.add(LSTM(64, input_shape = (PREVISOES, 1), return_sequences = False))
                model.add(Dense(1, kernel_initializer = initializer))
                model.add(BatchNormalization())
                model.add(Dense(1, kernel_initializer = initializer))
                
                # Configurando o EarlyStopping para o modelo não treinar mais que 3x seguidas se não obtiver melhorias nos resultados
                early_stopping = EarlyStopping(monitor = 'loss', patience = 3, mode = 'min')
                
                # Treinando o modelo com otimizador Adam
                model.compile(loss = 'mse', optimizer = Adam(learning_rate = learning_rate), metrics = ['mae'])
                h = model.fit(train_gen, epochs = EPOCAS, callbacks = [early_stopping], verbose = False)
                pred = np.empty(test.shape[0])
                
                # Realizando predições no conjunto de teste
                for i in range(len(pred)):
                    prediction = model.predict(test_gen, verbose = False)
                    pred[i] = prediction
                    test_input = np.delete(test_input, 0, 0)
                    test_input = np.concatenate((test_input, np.array(prediction).reshape(-1, 1)), axis = 0)
                    test_gen = TimeseriesGenerator(test_input, test_input, length = PREVISOES, batch_size = BATCH_SIZE)
                
                # Reorganizando o shape e chamando a função de calculo
                pred = ss.inverse_transform(pred.reshape(-1,1)).reshape(-1)
                pred, test, _ = calcFunction(train, pred, test)      

                mse.append(mean_squared_error(pred, test))
                
        except:
            mse.append(-1)
        
    return np.mean(mse)


# Report do algoritmo LSTM
def reportLSTM(data, modelName, calcFunction, learning_rate):
    global result
    global figs
    
    mse = []
    rmse = []
    mae = []
    mrv = []
    mpv = []
    
    title = modelName + ' - Long Short Term Memory (LSTM)'
    
    # Crossvalidation para cada parte do conjunto
    for train_id, test_id in split_range.split(data):
        train, test = data.iloc[train_id], data.iloc[test_id]
    
        # Normalização e reshape do conjunto de treino
        ss = StandardScaler()
        ss.fit(train.values.reshape(-1, 1))
        train_input = ss.transform(train.values.reshape(-1, 1))

        # Gerando conjunto de treino com TimeseriesGenerator baseado no conjunto atual e o batch informado
        test_input = train_input[-(PREVISOES + 1):]
        test_gen = TimeseriesGenerator(test_input, test_input, length = PREVISOES, batch_size = BATCH_SIZE)
        train_gen = TimeseriesGenerator(train_input, train_input, length = PREVISOES, batch_size = BATCH_SIZE)

        # Reset da sessão
        clear_session()
        
        # Construindo o modelo de LSTM com GlorotUniform pois inicializa de forma normalizada
        initializer = GlorotUniform(seed = SEED)
        model = Sequential()
        
        # 1 camada de LSTM com 64 entradas, 2 camadas densas e uma de normalização intermediando as camadas densas
        model.add(LSTM(64, input_shape = (PREVISOES, 1), return_sequences = False))
        model.add(Dense(1, kernel_initializer = initializer))
        model.add(BatchNormalization())
        model.add(Dense(1, kernel_initializer = initializer))
        
        # Configurando o EarlyStopping para o modelo não treinar mais que 3x seguidas se não obtiver melhorias nos resultados
        early_stopping = EarlyStopping(monitor = 'loss', patience = 3, mode = 'min')
        
        # Treinando o modelo com otimizador Adam
        model.compile(loss = 'mse', optimizer = Adam(learning_rate = learning_rate), metrics = ['mae'])
        h = model.fit(train_gen, epochs = EPOCAS, callbacks = [early_stopping], verbose = False)
        pred = np.empty(test.shape[0])

        # Realizando predições no conjunto de teste
        for i in range(len(pred)):
            prediction = model.predict(test_gen, verbose = False)
            pred[i] = prediction
            test_input = np.delete(test_input, 0, 0)
            test_input = np.concatenate((test_input, np.array(prediction).reshape(-1, 1)), axis = 0)
            test_gen = TimeseriesGenerator(test_input, test_input, length = PREVISOES, batch_size = BATCH_SIZE)

        # Reorganizando o shape e chamando a função de calculo
        pred = ss.inverse_transform(pred.reshape(-1,1)).reshape(-1)
        pred, test, indexPlot = calcFunction(train, pred, test)

        mse.append(mean_squared_error(pred, test, squared = True))
        rmse.append(mean_squared_error(pred, test, squared = False))
        mae.append(mean_absolute_error(pred, test))
        mrv.append(np.mean(test))
        mpv.append(np.mean(pred))
    
    result = record(result, title, np.mean(mse), np.mean(rmse), np.mean(mae), np.mean(mrv), np.mean(mpv), False)
    return plot(test.index, pred, mse, title, figs, axs[indexPlot], modelName)


# Gerando tabela de hypertunning com taxas de learning_rate
learning_rates = np.logspace(-5, 1, 7)
tuning = pd.DataFrame(learning_rates, columns=['learning_rate'])
tuning.head()


%%time
tuning['mse'] = tuning.apply(lambda x: GSLSTM(seriesHistory.copy(), x.learning_rate, calcPredTestOriginal), axis=1)

Wall time: 3min 15s


tuning.query('mse == mse.min() and mse != -1')


# Reset da figura
figs, axs = plt.subplots(nrows=4, sharex=True, figsize=(13,6))
figs.align_ylabels()
figs.tight_layout()
plt.close()


reportLSTM(seriesHistory.copy(), 'Original', calcPredTestOriginal, 0.0001)


learning_rates = np.logspace(-5, 1, 7)
tuning = pd.DataFrame(learning_rates, columns=['learning_rate'])
tuning.head()


%%time
tuning['mse'] = tuning.apply(lambda x: GSLSTM(deseasonal.copy(), x.learning_rate, calcPredTestDeseasonal), axis=1)

Wall time: 3min 18s


tuning.query('mse == mse.min() and mse != -1')


reportLSTM(deseasonal.copy(), 'Deseasonal', calcPredTestDeseasonal, 0.0001)


learning_rates = np.logspace(-5, 1, 7)
tuning = pd.DataFrame(learning_rates, columns=['learning_rate'])
tuning.head()


%%time
tuning['mse'] = tuning.apply(lambda x: GSLSTM(bc_history.copy(), x.learning_rate, calcPredTestBoxCox), axis=1)

Wall time: 3min 21s


tuning.query('mse == mse.min() and mse != -1')


reportLSTM(bc_history.copy(), 'BoxCox', calcPredTestBoxCox, 0.0001)


learning_rates = np.logspace(-5, 1, 7)
tuning = pd.DataFrame(learning_rates, columns=['learning_rate'])
tuning.head()


%%time
tuning['mse'] = tuning.apply(lambda x: GSLSTM(diff_history.copy(), x.learning_rate, calcPredTestStationary), axis=1)

Wall time: 3min 15s


tuning.query('mse == mse.min() and mse != -1')


reportLSTM(diff_history.copy(), 'Stationary', calcPredTestStationary, 0.0001)


result


# Tratando nomes e criando colunas de controle para os resultados gerados
topResult = (
    result 
    .assign(Full_Name = lambda x: x.Algorithm.apply(lambda a: a.split('(')[0]
                                                   .replace('ARIMA', 'Auto Arima')
                                                   .replace('Long Short Term Memory', 'LSTM')))
    .assign(Data_Category = lambda x: x.Algorithm.apply(lambda a: a.split(' - ')[0]))
    .assign(Algorithm = lambda x: x.Algorithm.apply(lambda a: a.split(' - ')[1].split('(')[0]
                                                   .replace('ARIMA', 'Auto Arima')
                                                   .replace('Long Short Term Memory', 'LSTM')))
    .sort_values('MSE')
)


topResult.head()


# Plot dos resultados obtidos por ordem ascendente do MSE

colors = {'Time Series Regression':'red',
          'Exponential Smoothing':'orange',
          'Auto Arima': 'green',
          'LSTM ': 'blue'}

# plotly figure
fig = go.Figure(layout = go.Layout(yaxis = {"type": "category", 'title': 'Algoritmo e Categoria'},
                        xaxis = {'title': 'MSE'},
                        title = 'MSE por Algoritmo e Tipo de Dado'))

for t in topResult['Algorithm'].unique():
    topResultFiltered = topResult[topResult['Algorithm']==t].copy()
    fig.add_traces(go.Bar(x = topResultFiltered['MSE'], y = topResultFiltered['Full_Name'], name = str(t),\
                          marker_color = str(colors[t]), orientation = 'h',
                          text = round(topResultFiltered['MSE']), textposition = 'outside'))
    
    
fig.update_layout(yaxis={'categoryorder':'total descending'}, autosize = False,
                  width = 1450,
                  height = 800)    
    
fig.show()

	index	order_purchase_timestamp
0	2016-09-15	1
1	2016-09-16	0
2	2016-09-17	0
3	2016-09-18	0
4	2016-09-19	0
...	...	...
709	2018-08-25	69
710	2018-08-26	73
711	2018-08-27	66
712	2018-08-28	39
713	2018-08-29	11

	Algorithm	MSE	RMSE	MAE	Mean_Real_Value	Mean_Predict_Value
0	Original - Time Series Regression	5485.688752	63.032946	45.057261	179.560484	164.066513
1	Deseasonal - Time Series Regression	5070.936469	57.842363	39.592783	179.560484	166.267257
2	BoxCox - Time Series Regression	5531.781002	63.453289	45.576572	179.560484	160.984010
3	Stationary - Time Series Regression	20765.011179	123.412496	88.116806	181.296371	175.123977

	Algorithm	MSE	RMSE	MAE	Mean_Real_Value	Mean_Predict_Value
0	Original - Time Series Regression	5485.688752	63.032946	45.057261	179.560484	164.066513
1	Deseasonal - Time Series Regression	5070.936469	57.842363	39.592783	179.560484	166.267257
2	BoxCox - Time Series Regression	5531.781002	63.453289	45.576572	179.560484	160.984010
3	Stationary - Time Series Regression	20765.011179	123.412496	88.116806	181.296371	175.123977
4	Original - Exponential Smoothing	5523.286110	60.505490	41.256041	179.560484	170.716575
5	Deseasonal - Exponential Smoothing	4995.372698	57.454750	38.814970	179.560484	167.407451
6	Original - Exponential Smoothing	5523.286110	60.505490	41.256041	179.560484	170.716575
7	Deseasonal - Exponential Smoothing	4995.372698	57.454750	38.814970	179.560484	167.407451
8	BoxCox - Exponential Smoothing	5295.550187	60.651937	41.800642	179.560484	174.995807
9	Stationary - Exponential Smoothing	48372.797213	196.181758	141.789213	181.296371	178.154334

	Algorithm	MSE	RMSE	MAE	Mean_Real_Value	Mean_Predict_Value
0	Original - Time Series Regression	5485.688752	63.032946	45.057261	179.560484	164.066513
1	Deseasonal - Time Series Regression	5070.936469	57.842363	39.592783	179.560484	166.267257
2	BoxCox - Time Series Regression	5531.781002	63.453289	45.576572	179.560484	160.984010
3	Stationary - Time Series Regression	20765.011179	123.412496	88.116806	181.296371	175.123977
4	Original - Exponential Smoothing	5523.286110	60.505490	41.256041	179.560484	170.716575
5	Deseasonal - Exponential Smoothing	4995.372698	57.454750	38.814970	179.560484	167.407451
6	BoxCox - Exponential Smoothing	5295.550187	60.651937	41.800642	179.560484	174.995807
7	Stationary - Exponential Smoothing	48372.797532	196.181759	141.789213	181.296371	178.154333
8	Original - ARIMA(2,1,1)(1,0,1)[5]	5734.428304	64.119816	46.036502	179.560484	168.922406
9	Deseasonal - ARIMA(2,1,1)(0,0,0)[0] intercept	9591.508344	77.052176	59.641647	179.560484	191.565756
10	BoxCox - ARIMA(4,1,2)(0,0,0)[5] intercept	11293.369481	88.024735	69.332624	179.560484	190.689776
11	Stationary - ARIMA(3,0,3)(0,0,0)[0] intercept	21477.900281	124.896247	86.287046	181.296371	160.562751

	learning_rate
0	0.00001
1	0.00010
2	0.00100
3	0.01000
4	0.10000

1. Preparando Conjunto de Pedidos¶

2. Iniciando Analise Seasonal¶

2.1 Testes¶

2.1.1 ADF:¶

2.1.2 KPSS:¶

3. Modelagem¶

3.1 TSR (Time Series Regression)¶

3.2 STL Deseasonal - TSR (Time Series Regression)¶

3.3 BoxCox - TSR (Time Series Regression)¶

3.4 Stationary - TSR (Time Series Regression)¶

3.5 Exponential Smoothing¶

3.6 Deseasonal - Exponential Smoothing¶

O código abaixo é uma replicação do item 3.5, de forma que só foi alterado a base de entrada de seriesHistory para deseasonal, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.7 BoxCox - Exponential Smoothing¶

O código abaixo é uma replicação do item 3.5, de forma que só foi alterado a base de entrada de seriesHistory para boxcox, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.8 Seasonal - Exponential Smoothing¶

O código abaixo é uma replicação do item 3.5, de forma que só foi alterado a base de entrada de seriesHistory para diff_history, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.9 ARIMA¶

3.10 Deseasonal - ARIMA¶

O código abaixo é uma replicação do item 3.9, de forma que só foi alterado a base de entrada de seriesHistory para deseasonal, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.11 BoxCox - ARIMA¶

O código abaixo é uma replicação do item 3.9, de forma que só foi alterado a base de entrada de seriesHistory para boxcox, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.12 Stationary - ARIMA¶

O código abaixo é uma replicação do item 3.9, de forma que só foi alterado a base de entrada original para stationary, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.13 LSTM¶

3.14 Deseasonal - LSTM¶

O código abaixo é uma replicação do item 3.13, de forma que só foi alterado a base de entrada original para deseasonal, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.15 BoxCox - LSTM¶

O código abaixo é uma replicação do item 3.13, de forma que só foi alterado a base de entrada original para boxcox, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

3.16 Stationary - LSTM¶

O código abaixo é uma replicação do item 3.13, de forma que só foi alterado a base de entrada original para stationary, assim verificando as diferenças de resultados ao utilizar diferentes transformações na base. Dessa forma, não terá comentários nesse item.¶

4. Comparação¶

	trend	season	periods	bias	method
0	add	add	5	True	L-BFGS-B
1	add	add	5	True	TNC
2	add	add	5	True	SLSQP
3	add	add	5	True	Powell
4	add	add	5	True	trust-constr

	Algorithm	MSE	RMSE	MAE	Mean_Real_Value	Mean_Predict_Value	Full_Name	Data_Category
5	Exponential Smoothing	4995.372698	57.454750	38.814970	179.560484	167.407451	Deseasonal - Exponential Smoothing	Deseasonal
1	Time Series Regression	5070.936469	57.842363	39.592783	179.560484	166.267257	Deseasonal - Time Series Regression	Deseasonal
13	LSTM	5138.516158	59.463829	41.561660	179.560484	160.851154	Deseasonal - LSTM	Deseasonal
15	LSTM	5141.851708	58.975721	40.576340	179.560484	178.591350	Stationary - LSTM	Stationary
6	Exponential Smoothing	5295.550187	60.651937	41.800642	179.560484	174.995807	BoxCox - Exponential Smoothing	BoxCox

	alpha	beta	gamma
0	0.0	0.0	0.0
1	0.0	0.0	0.1
2	0.0	0.0	0.2
3	0.0	0.0	0.3
4	0.0	0.0	0.4

	alpha	beta	gamma
0	0.0	0.0	0.0
1	0.0	0.0	0.1
2	0.0	0.0	0.2
3	0.0	0.0	0.3
4	0.0	0.0	0.4

	alpha	beta	gamma
0	0.0	0.0	0.0
1	0.0	0.0	0.1
2	0.0	0.0	0.2
3	0.0	0.0	0.3
4	0.0	0.0	0.4

	alpha	beta	gamma
0	0.0	0.0	0.0
1	0.0	0.0	0.1
2	0.0	0.0	0.2
3	0.0	0.0	0.3
4	0.0	0.0	0.4